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Abstract. We investigate finite particle systems of cold atoms bound in a local 
potential V(r). We derive the ground state energy and the particle density using 
a recently developed semiclassical theory (2008 Phys. Rev. Lett. 100 200408), 
and assuming the particles are described by the Haldane-Wu fractional exclusion 
statistics (FES) at unitarity. This approach is applied to atoms trapped into 
a three dimensional harmonic oscillator. We show that the parameter-free FES 
semiclassical theory yields results that are consistent with numerical simulations 
by Chang and Bertsch [2007 Phys. Rev A 76 021603(R)] and Bulgac (2007 Phys. 
Rev. A 76 040502). 



PACS numbers: 03.65.Sq, 03.75.Ss, 05.30.Fk, 71.10.-w 

Recently, the crossover between the BEC (Bosc-Einstcin condensate) and 
BCS (Bardcen-Cooper-SchricfTcr superconductivity) regimes in a cold atom gas has 
attracted a great deal of attention 1 . Using the effect of the Feshbach resonances, 
the interactions or the scattering length a can be tuned by an external magnetic 
field. When increasing the magnetic field, a grows from a small negative value in 
the BCS regime to a small positive value in the condensate state. In between these 
two regimes, the scattering length diverges and the gas is said to be at unitarity. 
Since no relevant energy scale exists but the Fermi wave length, one expects a 
universal behavior of the gas [2]. In that picture Papenbrock [3] noticed that the 
ground state energy of this system grows like the Thomas-Fermi (TF) approximation. 
The proportionality constant, expected to be universal, has been recently computed 
analytically by Bhaduri et al. [HOH], assuming interactions at unitarity are simulated 
by non-interacting particles which obey the fractional Haldane-Wu's statistics (FES) 
[7] . This extension of standard Bose-Einstein and Fermi-Dirac statistics was computed 
in connection with the fractional quantum Hall effect, and later extended by Wu [SJ. 
In FES, the Pauli exclusion principle is generalized such that the occupation number 
follows n g (e, A) = l/(w[e( f -~ x ^ kBT ]+g) where T is the temperature, g is the occupancy 
factor, A is the chemical potential and w 9 (x)[l + w{x)] 1 ~ 9 = x. We recover the Bosc 
and Fermi distributions, respectively, with g = and g = 1. At zero temperature the 
FES occupation number obeys the step distribution 

n g (e,X) = - 6(e-A) , (1) 
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where Q(x) is the Heaviside function. For such statistics the semiclassical theory, 
initiated by the trace formula of Gutzwillcr [5j and recently extended in |10j for 
spatially-varying densities, provides a convenient formalism to compute the relevant 
thermodynamic quantities for finite particle systems [TTJ EH Q3J . 

In this communication, we establish the FES semiclassical theory for the ground 
state energy and for the particle density. We discuss their expectation values for 
amplitudes and shell effects. We check to what extent the FES hypothesis of unitary 
atoms trapped into a three dimensional isotropic harmonic potential (IHO) agrees with 
benchmark numerical simulations performed by Chang and Bertsch j!5) and Bulgac 

m- 

We first derive the semiclassical theory for N non-interacting particles with 
mass m, obeying FES statistics and bound by a local smooth potential V(r) at 
zero temperature. Here we restrict ourselves to three-dimensional systems, but the 
generalization to other dimensions is straightforward. The discrete energy eigenvalues 
ej and cigenfunctions ipji?) are given by the stationary Schrodingcr equation. Using 
the FES occupation number ((T|), the quantum- mechanical single-particle (SP) level 
density and particle density of the system at zero temperature, are given respectively 
by K e ) = E £j S ( e ~ e i) and 

P (r) = 2^ % ( £j , A)V*(r)^-(r) = ~T J ^W^W . (2) 

The chemical potential is determined by the integrated level density which counts the 
number of SP state up to A: 

/■OO r\ i> A 

Af(X) = 2 n g (e,\)v(e)de= - / v{e)de. (3) 
Jo 9 Jo 

The factor 2 in ((2]) and (J3j) takes into account the spin degeneracy. For a given 

potential the SP level density and the particle density are decomposed into a smooth 

part v (resp. p) coming from the TF theory plus an oscillatory contribution 8v (resp. 

8p): 

v{e)=V{e) + 5v{e), (4) 
p(r) = p(r) + 6p(r). (5) 

Below, we will show in more detail how to describe these two components. (j4]) induces 
a similar decomposition for the integrated level density: 

Af(\) = Af(X) + 5M(X) 

2 f x 2 f x 

= - / V{e) de +- 6v(e) de . (6) 
9 Jo 9 Jo 

In order to study a system with a finite number of particles, we use the canonical 
expressions for thermodynamic quantities. For a system of N non-interacting 
fcrmions, we define its ground state energy E(N), the shell correction energy SE(N) 
and the smooth TF component E(N) as [TTllTg] : 

SE(N) = E(N) - E(N) 



oo 



2 / n g (e, \)ev(e)de — 2 / n g (e, \)ev(e)de 
Jo Jo 

2 f x 2 f x 

- / tv(e)dt - - / ev(e)de. (7) 
9 Jo 9 Jo 
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The chemical potential A and its smooth part A fix the number of particles. They are 
defined by inversion of the exact and the average integrated level densities, N(\) = N 
and 

Af(\) = N , (8) 

respectively. Note that, from ([6]) and ((HJ, both A and A depend on g. Unfortunately, (0 
is difficult to exploit analytically because the discretization of A is difficult to impose. 
From ([7]) it can be shown that, neglecting terms of second order in the parameter 
A — A, SE may be approximated by [19 l [20j: 



SE(N) w - / 5N{e)de . (9) 
Jo 

This, together with the definitions of E(N), p(r) and 5p(r), are the basic equations 
upon which we will base our analysis of the ground state of FES Fermi gases. 

In the limit N — > oo, p and E are expected to carry over into the approximation 
obtained in the TF theory [21]. v is given for any local potential V(r), by 

m 3/2 r 

v = v,* = ^-7= J y/e-V(q)e[e - y(q)]dq. (10) 

The TF expression for the particle density has to be weighted by the FES occupation 
number (pj following (|5J). It yields: 

l2 [\-V{r)f'\ (11) 



where A is given by inverting (j5J) with the l.h.s. of (jSJ) and (|10[) . E is computed using 
the r.h.s. of © and (JTDJ : 

= n2h3y /2 J J '\/e-V(<i)e[e Vi^dqde. (12) 

For the oscillating parts Sp and 5E, we now use the main formulas from [91 [TUI [T^] 
for the special case at D = 3 space dimensions. Including the FES occupancy factor, 
the semiclassical expression for dp, to leading order in h, yields 



^ cos -^-^--tt . (13) 



The sum is, in general, over all non-periodic orbits 7 starting and ending in r taken 
at energy A. 5 7 is the action function 5 7 = J r r _r p(A, q) • dq, where p(A,q) is the 
classical momentum in the point q. /i 7 is the Morse index that counts the number of 
conjugate points along the orbit [9]. X>j_ 7 = (9px/9r^)|r'=r is the stability determinant 
calculated from the components p± and transverse to the orbit 7 of the initial 
momentum and final coordinate, respectively. Here p is the modulus of the momentum 
in r and T 7 = dS' 7 (e, r)/de| e= ^ is the running time of the orbit 7. 

To leading order in h, the oscillating part 5v(e) is given by the semiclassical trace 
formula 



5v(e) ~ VS PO (e) cos 



PO 



(14) 



where the sum runs over all periodic orbits (POs). For systems in which all orbits 
are isolated in phase space, Gutzwiller [9] derived explicit expressions for the smooth 
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amplitudes B PO (e), which depend on the stability of the orbits, and for the Maslov 
indices a PO . Performing the trace integral in the semiclassical Green function [5] along 
all directions transverse to each orbit 7, the stationary phase approximation (SPA) 
leads immediately to the periodicity of the contributing orbits. The Maslov index er po 
collects all phases occurring in the semiclassical Green function and in the SPA for the 
trace integral (see [22] for detailed computations of cr PO ). S PO (e) is the closed action 
integral S po (e) = j> pQ p(e, q) • dq. We compute analytically 5E using the r.h.s. of ((6]) 
and (J9j> . To leading order in H we get 

2ft 2 ^£ po (A)_[l Cl ~, tt_ 1 



SE(N) = > pov J y cos 

9 ^ T 2 (A) 



h S PQ {\)--o P 



where T PO is the running time of the PO. 

We emphasize that the semiclassical approximations are not valid in regions close to 
the classical turning points r x defined by V(r\) = A. Since the classical momentum 
p becomes zero, the spatial density (| X3[) always diverges at the turning points. 
Furthermore the running time which appears in the denominator of (|13[) , may vanish at 
the turning point for certain orbits. These divergences can be overcome by linearizing 
the smooth potential V(r) around the classical turning points [231 l2~i] . Note that 
the phases and amplitudes in (fT3f and (fT5"j) strongly depend on g through A, but no 
simple behavior clearly emerges. In practice, one has to compute them explicitly for 
a given potential. In the following, we address this issue in the special case of a 3D 
IHO applied to a unitary Fermi gas. 

Bhaduri et al [H [5J [5] assume that FES simulates interactions between particles 
at unitarity and give the analytical expression of the unitary occupancy factor 
g u = 1 - V2/2 « 0.29 for a 3D IHO confinement V(r) = tocjV/2 with |r| = r. 
They compute observables in the extended TF limit leading to smooth variations 
only. Generalizing their method, we introduce oscillating corrections that depend on 
h for g — g u . For a 3D IHO, the smooth functions p and E are easily computed and 
give 

»" =3^(^) 3/2(X — V/2)3/2 ' <16) 



E(N) = ^ 



Ahuj) + 8 V Ttw 



■o(l) 



(17) 



with \{N) = hbj^g^N) 1 / 3 . These two results are consistent with those in [SJ [23] , 
Using the formulas in p~3l [26] for dp and 5E with the FES modifications ([T3]) and 
US]), we get: 



-ro 2 w(3JV) 1/a C0S { S ± 



- 2 ; 2 /3 2. v tW — - . 



1/3 2 7T 2 ^ k 2 

yu fc=i 



(19) 



Here we have used the analytical form of the actions and periods = 
(2k + l)ir\/u> =F rp =p 2A/w arctan(ma;r/p), = (2fc + l)ir/u> =p 2/w arctan(ma;r/p), 
/Uj^^ = 6fc + 1 and fj,_ = 6fc + 3. From ([T9| we note that the amplitude of the shell 
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correction energy for a unitary Fermi gas is larger (5E/E oc g u ) compare to the 
standard Fermi gas (g = 1). The approximate frequency of the shell fluctuations 
is given by the phase of the cosine function of the k = 1 term in the sum (TIT)]) , 
Thus, the closed-shell (resp. mid-shell) numbers, given by the values of N that 
minimize (resp. maximize) 5E(N), are well approximated by N cs = i 3 /(3g u ) [resp. 
N m s = (i + l/2) 3 /(3<7u)] with i e N. The presence of p in the denominator of (|18j) 
gives no simple behavior for dp in g n . Nevertheless close to r s=s 0, the ratio Sp/p grows 

— 1/2 

like <7 U , so that the oscillations are amplified by a factor 1.9. 
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Figure 1. Upper panel: ground state energy of a Fermi gas at unitarity as a 
function of the particle number N in a 3D IHO. The solid line corresponds to the 
scmiclassical FES ground state energy l|17fl + l|19j l. the dotted (resp. dashed) line 
corresponds to the numerical data 1151 (resp. 1161 ) (units h = m . = u = 1). Lower 
panel: shell correction energy for the same system. The solid line corresponds to 
the semiclassical FES theory JT9}. The dotted (resp. dashed) line is the numerical 
data 1151 (resp. 1161 ) subtracted by the smooth energy H17I I. 

When r is close to the classical turning point r\ = uj^ 1 {2\/m) 1 / 2 , the 
semiclassical approximation (1181) diverges. Following the regularisation method 
detailed in , the density profile becomes 

= 7^{AiWAi'(z) + 2z[Ai'(z)] 2 -2z 2 Ai 2 (z)}, (20) 

where p = 2(2muj/h 2 ) 1/3 (2mX) 1/6 , z = p (r - r\/2)/2, Ai(x) is the Airy function 
and Ai'(a:) its derivative. 

We now compare ours results to the ab initio Green function Monte Carlo method 
(GFMC) from Chang and Bertsch [15] and the supcrfluid local density approximation 
(SLDA) computed by Bulgac [16]. Figure [1] focuses on the ground state energy for 
N = 2-22 atoms. The upper panel shows a qualitative agreement between the three 
curves. In the lower panel only the oscillating component is plotted. While the 
shell effects are observed for the FES model (solid line), they arc not present in the 
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Figure 2. Particle density of a Fermi gas at unitarity as a function of r for 
N = 20 particles in a 3D IHO. The solid line corresponds to the semiclassical 
FES particle density 116t + 118t . the dotted (resp. dashed) line corresponds to the 
numerical data |15| (resp. |16| ) (units H = m = ui = 1). The dottcd-dashed line 
is the semiclassical improvement discussed in the text. 

SLDA model (dashed line). The reason is LDA theories correspond to (extended) TF 
approximations [27]. They only provide the average part of thermodynamic quantities, 
but cannot take into account shell effects. The GFMC data (dotted line) shows some 
irregularities but no clear oscillations. This is not surprising since the semiclassical 
theories are known to be more accurate for large particle numbers. Hence further 
numerical data are needed in order to check the presence of shell effects at unitarity. 
Note the odd-even oscillation of the energy in the numerical computations, which is 
due to the pairing correlations not included in the FES theory. 

In figure El we plot the particle density for N = 20. The GFMC (dotted line) 
and FES (solid line) models show some density oscillations. Although the order of 
magnitude is correct, the oscillations are not in phase. This is attributed to the 
discrepancy of the semiclassical theory for particle numbers far from the shell or mid- 
shell closure, thus leading to an inaccurate sign of dp (see [24] for an exhaustive 
discussion). The FES model is expected to give best results for N « N cs and TV sa iV ms . 
The SLDA model (dashed line) gives an accurate average profile of the particle density 
but shows no oscillations. We mention that same comparisons have been done with 
the numerical work of von Stecher et al. (28[ [29] [30] [31] leading to similar results 
with the SLDA model. The linear behavior of the FES particle density observed for 
r > 1.5, is a consequence of the breakdown of the semiclassical approximation near 
the classical turning point. We recover the usual exponential tail using the regularised 
formula l|2"0|) . The dotted-dashed line illustrates an improved density profile for which 
we have changed the overall sign into (| 18[) and switched to the regularised expression 
(|20p for r > 1.5: in this case, the agreement with GFMC is much better. 

In conclusion, we have investigated finite fermions systems which are described 
by the Haldane-Wu statistics going beyond the TF approximation. The ground 
state energy and the particle density of this system are derived analytically at zero 
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temperature. We used the FES semiclassical theory as a parameter-free model of 
unitary Fermi gases and we discussed shell effects as a function of the particle number 
and the position. Considering that the semiclassical model is expected to give 
better results for large N, we gave reasonably good agreement with two numerical 
studies. The investigation of the more general occupation number distribution 
n g (e, A) = l/(w[e^ _A ^ fcBT ] + g) in the case of non-zero temperatures is let for future 
works. 

I acknowledge A Bulgac, S Y Chang and D Blume for providing the numerical 
data and M Brack for fruitful discussions and constant encouragement. This work 
was funded in part by the French National Research Agency ANR (project ANR-06- 
BLAN-0059). 
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